
include "stdsage.pxi"
include "cdefs.pxi"

cdef extern from "cmath":
    double floor(double)

#cdef extern from "stdlib.h":
#    void free(void* ptr)
#    void* malloc(size_t size)

from sage.rings.integer cimport Integer
from sage.rings.rational cimport Rational
from sage.structure.element cimport Element, FieldElement, RingElement, ModuleElement
from sage.structure.parent_base cimport ParentWithBase
from sage.libs.ntl.ntl_ZZX cimport ntl_ZZX
from sage.libs.ntl.ntl_ZZ cimport ntl_ZZ

from sage.rings.integer_ring import IntegerRing
from sage.rings.number_field.number_field import NumberField
from sage.schemes.elliptic_curves.constructor import EllipticCurve_from_c4c6
from psage.ellcurve.lseries.aplist_sqrt5 import aplist as compute_aplist

R = IntegerRing()['x']
f = R([-1, -1, 1])
K = NumberField(f, 'a')

from sage.libs.gmp.mpz cimport (mpz_t, mpz_set, mpz_set_si, mpz_init, mpz_clear, mpz_fdiv_ui)

from sage.rings.number_field.number_field_element cimport NumberFieldElement, NumberFieldElement_absolute

from sage.rings.number_field.number_field_element_quadratic cimport NumberFieldElement_quadratic

cdef NumberFieldElement_quadratic a = K.gen()

cdef looks_square(double x1, double x2):
    # checks if x looks like a square in K, where x1 and x2 are two embeddings of x
    #
    # assumes that it is an integer. this function might return true when the answer
    # is false, because i am being lazy, but it won't fail the other way around.

    if x1 < 0 or x2 < 0:
        return -1

    cdef double z1, z2, w1
    
    z1 = sqrt(x1)
    z2 = -z1
    w1 = sqrt(x2)

    cdef double trace1, trace2

    trace1 = z1 + w1
    trace2 = (z1 - w1)/sqrt(5.0)

    if abs(trace1 - floor(trace1 + .5)) < .0001 and abs(trace2 - floor(trace2 + .5)) < .0001:
        return 1

    trace1 = z2 + w1
    trace2 = (z2 - w1)/sqrt(5.0)
    if abs(trace1 - floor(trace1 + .5)) < .0001 and abs(trace2 - floor(trace2 + .5)) < .0001:
        return 1

    return 0

def search_possibilities(D_guesses,
                         D_guess_embeddings,
                         _j1a_possibilities,
                         _j1b_possibilities,
                         _j2a_possibilities,
                         _j2b_possibilities,
                         level,
                         aplist,
                         int verbose):

    # start by converting everything to C arrays, since Cython doesn't seem to do
    # this for us (not that it could for the first entry)

    cdef int num_D_guesses = len(D_guess_embeddings)
    #cdef double * D_embeddings1 = <double *>malloc(num_D_guesses * sizeof(double))
    #cdef double * D_embeddings2 = <double *>malloc(num_D_guesses * sizeof(double))
    cdef int n
    cdef int m
    #for n in range(num_D_guesses):
    #    D_embeddings1[n] = _D_guess_embeddings[n][0]
    #    D_embeddings2[n] = _D_guess_embeddings[n][1]
        
    cdef int num_j_guesses = len(_j1a_possibilities)
    cdef double * j1a_possibilities = <double *>malloc(num_j_guesses * sizeof(double))
    cdef double * j1b_possibilities = <double *>malloc(num_j_guesses * sizeof(double))
    cdef double * j2a_possibilities = <double *>malloc(num_j_guesses * sizeof(double))
    cdef double * j2b_possibilities = <double *>malloc(num_j_guesses * sizeof(double))
    for n in range(num_j_guesses):
        j1a_possibilities[n] = _j1a_possibilities[n]
        j1b_possibilities[n] = _j1b_possibilities[n]
        j2a_possibilities[n] = _j2a_possibilities[n]
        j2b_possibilities[n] = _j2b_possibilities[n]
    cdef NumberFieldElement_quadratic D, c4, c6, Z
    cdef double D1
    cdef double D2
    cdef double j1
    cdef double j2
    cdef double Z1
    cdef double Z2

    cdef double c4a, c4b

    cdef double * j1_possibilities
    cdef double * j2_possibilities
    cdef int count = 0
    
    cdef int total_guesses = num_D_guesses * num_j_guesses

    cdef int A
    cdef int B

    cdef int x
    cdef int y

    cdef int count2 = 0

    cdef int foo

    for n in range(num_D_guesses):
        D = D_guesses[n]
        D1 = D_guess_embeddings[n][0]
        D2 = D_guess_embeddings[n][1]

        if D1 > 0:
            j1_possibilities = j1a_possibilities
        else:
            j1_possibilities = j1b_possibilities
        if D2 > 0:
            j2_possibilities = j2a_possibilities
        else:
            j2_possibilities = j2b_possibilities 

        for n in range(num_j_guesses):
            j1 = j1_possibilities[n]
            j2 = j2_possibilities[n]
            if verbose:
                count += 1
                if count % 1000 == 0:
                    print "Completed", count, "guesses out of", total_guesses
                    print c4a, c4b, A, B
            if j1 * D1 > 0:
                c4a = (j1 * D1)**(1.0/3.0)
            else:
                c4a = -(-j1 * D1)**(1.0/3.0)
            if j2 * D2 > 0:
                c4b = (j2 * D2)**(1.0/3.0)
            else:
                c4b = -(-j2 * D2)**(1.0/3.0)

            #Z1 = c4a**3 - 1728 * D1
            #Z2 = c4b**3 - 1728 * D2
            
            A = <int>floor(c4a + c4b + .5)
            B = <int>floor( (c4a - c4b)/sqrt(5.0) + .5)
            
            #note that A and B are integers, and they ought to have the same
            #parity. if they don't then we just add one to B

            if A % 2 != B % 2:
                B += 1

            # our choice of c4 is going to be
            # A/2 + (B/2)\sqrt 5
            
            # now that we've rounded, we redefine c4a and c4b to be the embeddings of this choice
            # of c4

            c4a = A/2.0 + B/2.0 * sqrt(5.0)
            c4b = A/2.0 - B/2.0 * sqrt(5.0)

            if not looks_square(c4a**3 - 1728*D1, c4b**3 - 1728*D2):
                continue

            # this is where we construct number field elements for the first time
            c4 = a._new()
            if A % 2 == 1:
                mpz_set_si(c4.a, A)
                mpz_set_si(c4.b, B)
                mpz_set_si(c4.denom, 2)
            else:
                mpz_set_si(c4.a, A/2)
                mpz_set_si(c4.b, B/2)
                mpz_set_si(c4.denom, 1)



            Z = c4**3 - 1728 * D

            # this is where things can get slow
            for c6 in Z.sqrt(True):
                print "Trying c4 =", c4, "; c6 = ", c6, " ; D =", D
                E = EllipticCurve_from_c4c6(c4, c6).integral_model().global_minimal_model()
                if E.conductor() == level:
                    new_aplist = compute_aplist(E, 200)
                    if new_aplist == aplist[:len(new_aplist)]:
                        return E

            #if Z.is_square():
            #    print "square!"

    print count2, "false positives."
    return 1
